rm(list = ls())
gc()

library(lubridate)
library(tidyverse)
library(foreign)
library(ggplot2)
library(ggthemes)
library(estimatr)
library(mice)
library(ggpubr)
library(texreg)
df_long = readRDS("data/baseline_clean_mi.rds")

df = as.mids(df_long)

indices = df$data %>% 
  select(c(response_num, names(df$data)[str_detect(names(df$data), "^index")]))

covars = c("syr_origin", "syr_origin_urban", "location_its", "sick", "head_educ2", "toddler", "elderly", "female_headed_hh", "location_district", "hezb_control")
covars2 = c("syr_origin_urban", "location_its", "sick", "head_educ2", "toddler", "elderly", "female_headed_hh", "hezb_control")

source("functions/plot_it.R")
source("functions/effect_plotter3.R")

# Return ever--------------------------------------------

## Individual indices--------------------------------------------
predictors = names(indices)[2:length(indices)]

for(i in 1:length(predictors)){
  if(predictors[i] != "index_mobility"){
    assign(paste0("reg", i), with(df, lm_robust(as.formula(paste("rtrn_head_ever", paste(c(predictors[i], covars), collapse = " + "), sep = " ~ ")), clusters = location_town, weights = weights)) %>% pool() %>% summary(conf.int = TRUE) %>% mutate(outcome = "expect2yr_syria"))
  }else{
    assign(paste0("reg", i), with(df, lm_robust(as.formula(paste("rtrn_head_ever", paste(c(predictors[i], covars2), collapse = " + "), sep = " ~ ")), clusters = location_town, weights = weights)) %>% pool() %>% summary(conf.int = TRUE) %>% mutate(outcome = "expect2yr_syria"))
  }
}

dfplot1 = reg1

for(i in 2:length(predictors)){
  dfplot1 = dfplot1 %>% bind_rows(eval(parse(text = paste0("reg", i))))
}


dfplot1 = dfplot1 %>% mutate(bold = "plain") %>% 
  mutate(term = factor(term, levels = c("index_safety_syria", "index_regime_control",
                                        "index_econ_syria", "index_services_syria", "index_networks_syria",
                                        "index_econ_leb",
                                        "index_services_leb", "index_networks_lebanon", "index_soc_leb",
                                        "index_legal_leb", "index_mobility", "index_family",
                                        "index_info_quality"))) %>% 
  mutate(term = recode(term, "index_econ_leb" = "Economic well-being (Leb.)",
                       "index_soc_leb" = "Social well-being (Leb.)",
                       "index_services_leb" = "Services (Leb.)",
                       "index_legal_leb" = "Legal conditions (Leb.)",
                       "index_safety_syria" = "Safety (Syr.)",
                       "index_econ_syria" = "Economic well-being (Syr.)",
                       "index_services_syria" = "Services (Syr.)",
                       "index_networks_syria" = "Networks (Syr.)",
                       "index_mobility" = "Log travel distance", 
                       "index_family" = "Log household size",
                       "index_networks_lebanon" = "Networks (Leb.) ",
                       "index_info_quality" = "Confidence in information",
                       "index_regime_control" = "Regime control (Syr.)")) %>% 
  #mutate(term = fct_rev(term)) %>% 
  mutate(outcome = recode(outcome, "rtrn_head_12" = "Return in 12 months",
                          "rtrn_head_ever" = "Return ever",
                          "prep_pca" = "Prepare to return", 
                          "expect2yr_syria" = "Return in 2 years"),
         outcome = factor(outcome, levels = c("Return in 12 months", "Prepare to return", "Return in 2 years", 
                                              "Return ever"))
  )

dfplot1 = dfplot1 %>% arrange(outcome, term)
dfplot1 = dfplot1 %>% filter(!is.na(term))


## All indices--------------------------------------------
predictors = names(indices)[2:length(indices)]
predictors_a = predictors[which(predictors!="index_services_syria")]
predictors_b = predictors[which(predictors!="index_safety_syria")]

reg1a = with(df, lm_robust(as.formula(paste("rtrn_head_ever", paste(c(predictors_a, covars2), collapse = " + "), sep = " ~ ")), weights = weights, clusters = location_town)) %>% pool %>% summary(conf.int = TRUE) %>% mutate(outcome = "expect2yr_syria")
reg1b = with(df, lm_robust(as.formula(paste("rtrn_head_ever", paste(c(predictors_b, covars2), collapse = " + "), sep = " ~ ")), weights = weights, clusters = location_town)) %>% pool %>% summary(conf.int = TRUE) %>% mutate(outcome = "expect2yr_syria")
reg1 = reg1a %>% bind_rows(reg1b %>% filter(term == "index_services_syria"))

dfplot2 = reg1

dfplot2 = dfplot2 %>% filter(term %in% names(indices))

dfplot2 = dfplot2 %>% mutate(bold = "plain") %>% 
  mutate(term = factor(term, levels = c("index_safety_syria", "index_regime_control",
                                        "index_econ_syria", "index_services_syria", "index_networks_syria",
                                        "index_econ_leb",
                                        "index_services_leb", "index_networks_lebanon", "index_soc_leb",
                                        "index_legal_leb", "index_mobility", "index_family",
                                        "index_info_quality"))) %>% 
  mutate(term = recode(term, "index_econ_leb" = "Economic well-being (Leb.)",
                       "index_soc_leb" = "Social well-being (Leb.)",
                       "index_services_leb" = "Services (Leb.)",
                       "index_legal_leb" = "Legal conditions (Leb.)",
                       "index_safety_syria" = "Safety (Syr.)",
                       "index_econ_syria" = "Economic well-being (Syr.)",
                       "index_services_syria" = "Services (Syr.)",
                       "index_networks_syria" = "Networks (Syr.)",
                       "index_mobility" = "Log travel distance", 
                       "index_family" = "Log household size",
                       "index_networks_lebanon" = "Networks (Leb.) ",
                       "index_info_quality" = "Confidence in information",
                       "index_regime_control" = "Regime control (Syr.)")) %>% 
  #mutate(term = fct_rev(term)) %>% 
  mutate(outcome = recode(outcome, "rtrn_head_12" = "Return in 12 months",
                          "rtrn_head_ever" = "Return ever",
                          "prep_pca" = "Prepare to return", 
                          "expect2yr_syria" = "Return in 2 years"),
         outcome = factor(outcome, levels = c("Return in 12 months", "Prepare to return", "Return in 2 years", 
                                              "Return ever"))
  )


dfplot2 = dfplot2 %>% arrange(outcome, term)

mod1 = createTexreg(coef.names = as.character(dfplot1$term), coef = dfplot1$estimate, se = dfplot1$std.error, pvalue = dfplot1$p.value, model.name = "Ever (Ind. Indices)")
mod2 = createTexreg(coef.names = as.character(dfplot2$term), coef = dfplot2$estimate, se = dfplot2$std.error, pvalue = dfplot2$p.value, model.name = "Ever (All Indices)")



# Return in 2 years--------------------------------------------
## Individual indices--------------------------------------------
predictors = names(indices)[2:length(indices)]

for(i in 1:length(predictors)){
  if(predictors[i] != "index_mobility"){
    assign(paste0("reg", i), with(df, lm_robust(as.formula(paste("expect2yr_syria", paste(c(predictors[i], covars), collapse = " + "), sep = " ~ ")), clusters = location_town, weights = weights)) %>% pool() %>% summary(conf.int = TRUE) %>% mutate(outcome = "expect2yr_syria"))
  }else{
    assign(paste0("reg", i), with(df, lm_robust(as.formula(paste("expect2yr_syria", paste(c(predictors[i], covars2), collapse = " + "), sep = " ~ ")), clusters = location_town, weights = weights)) %>% pool() %>% summary(conf.int = TRUE) %>% mutate(outcome = "expect2yr_syria"))
  }
}

dfplot1 = reg1

for(i in 2:length(predictors)){
  dfplot1 = dfplot1 %>% bind_rows(eval(parse(text = paste0("reg", i))))
}


dfplot1 = dfplot1 %>% mutate(bold = "plain") %>% 
  mutate(term = factor(term, levels = c("index_safety_syria", "index_regime_control",
                                        "index_econ_syria", "index_services_syria", "index_networks_syria",
                                        "index_econ_leb",
                                        "index_services_leb", "index_networks_lebanon", "index_soc_leb",
                                        "index_legal_leb", "index_mobility", "index_family",
                                        "index_info_quality"))) %>% 
  mutate(term = recode(term, "index_econ_leb" = "Economic well-being (Leb.)",
                       "index_soc_leb" = "Social well-being (Leb.)",
                       "index_services_leb" = "Services (Leb.)",
                       "index_legal_leb" = "Legal conditions (Leb.)",
                       "index_safety_syria" = "Safety (Syr.)",
                       "index_econ_syria" = "Economic well-being (Syr.)",
                       "index_services_syria" = "Services (Syr.)",
                       "index_networks_syria" = "Networks (Syr.)",
                       "index_mobility" = "Log travel distance", 
                       "index_family" = "Log household size",
                       "index_networks_lebanon" = "Networks (Leb.) ",
                       "index_info_quality" = "Confidence in information",
                       "index_regime_control" = "Regime control (Syr.)")) %>% 
  #mutate(term = fct_rev(term)) %>% 
  mutate(outcome = recode(outcome, "rtrn_head_12" = "Return in 12 months",
                          "rtrn_head_ever" = "Return ever",
                          "prep_pca" = "Prepare to return", 
                          "expect2yr_syria" = "Return in 2 years"),
         outcome = factor(outcome, levels = c("Return in 12 months", "Prepare to return", "Return in 2 years", 
                                              "Return ever"))
  )

dfplot1 = dfplot1 %>% arrange(outcome, term)
dfplot1 = dfplot1 %>% filter(!is.na(term))


## All indices--------------------------------------------
predictors = names(indices)[2:length(indices)]
predictors_a = predictors[which(predictors!="index_services_syria")]
predictors_b = predictors[which(predictors!="index_safety_syria")]

reg1a = with(df, lm_robust(as.formula(paste("expect2yr_syria", paste(c(predictors_a, covars2), collapse = " + "), sep = " ~ ")), weights = weights, clusters = location_town)) %>% pool %>% summary(conf.int = TRUE) %>% mutate(outcome = "expect2yr_syria")
reg1b = with(df, lm_robust(as.formula(paste("expect2yr_syria", paste(c(predictors_b, covars2), collapse = " + "), sep = " ~ ")), weights = weights, clusters = location_town)) %>% pool %>% summary(conf.int = TRUE) %>% mutate(outcome = "expect2yr_syria")
reg1 = reg1a %>% bind_rows(reg1b %>% filter(term == "index_services_syria"))

dfplot2 = reg1

dfplot2 = dfplot2 %>% filter(term %in% names(indices))

dfplot2 = dfplot2 %>% mutate(bold = "plain") %>% 
  mutate(term = factor(term, levels = c("index_safety_syria", "index_regime_control",
                                        "index_econ_syria", "index_services_syria", "index_networks_syria",
                                        "index_econ_leb",
                                        "index_services_leb", "index_networks_lebanon", "index_soc_leb",
                                        "index_legal_leb", "index_mobility", "index_family",
                                        "index_info_quality"))) %>% 
  mutate(term = recode(term, "index_econ_leb" = "Economic well-being (Leb.)",
                       "index_soc_leb" = "Social well-being (Leb.)",
                       "index_services_leb" = "Services (Leb.)",
                       "index_legal_leb" = "Legal conditions (Leb.)",
                       "index_safety_syria" = "Safety (Syr.)",
                       "index_econ_syria" = "Economic well-being (Syr.)",
                       "index_services_syria" = "Services (Syr.)",
                       "index_networks_syria" = "Networks (Syr.)",
                       "index_mobility" = "Log travel distance", 
                       "index_family" = "Log household size",
                       "index_networks_lebanon" = "Networks (Leb.) ",
                       "index_info_quality" = "Confidence in information",
                       "index_regime_control" = "Regime control (Syr.)")) %>% 
  #mutate(term = fct_rev(term)) %>% 
  mutate(outcome = recode(outcome, "rtrn_head_12" = "Return in 12 months",
                          "rtrn_head_ever" = "Return ever",
                          "prep_pca" = "Prepare to return", 
                          "expect2yr_syria" = "Return in 2 years"),
         outcome = factor(outcome, levels = c("Return in 12 months", "Prepare to return", "Return in 2 years", 
                                              "Return ever"))
  )


dfplot2 = dfplot2 %>% arrange(outcome, term)



mod3 = createTexreg(coef.names = as.character(dfplot1$term), coef = dfplot1$estimate, se = dfplot1$std.error, pvalue = dfplot1$p.value, model.name = "2 years (Ind. Indices)")
mod4 = createTexreg(coef.names = as.character(dfplot2$term), coef = dfplot2$estimate, se = dfplot2$std.error, pvalue = dfplot2$p.value, model.name = "2 years (All Indices)")


# Robustness check using the services regressions ----------------
predictors = names(indices)[2:length(indices)]
predictors_a = predictors[which(predictors!="index_services_syria")]
predictors_b = predictors[which(predictors!="index_safety_syria")]

reg1a = with(df, lm_robust(as.formula(paste("rtrn_head_12", paste(c(predictors_a, covars2), collapse = " + "), sep = " ~ ")), weights = weights, clusters = location_town)) %>% pool %>% summary(conf.int = TRUE) %>% mutate(outcome = "rtrn_head_12")
reg1b = with(df, lm_robust(as.formula(paste("rtrn_head_12", paste(c(predictors_b, covars2), collapse = " + "), sep = " ~ ")), weights = weights, clusters = location_town)) %>% pool %>% summary(conf.int = TRUE) %>% mutate(outcome = "rtrn_head_12")
reg1 = reg1b %>% bind_rows(reg1a %>% filter(term == "index_safety_syria"))

reg2a = with(df, lm_robust(as.formula(paste("rtrn_head_ever", paste(c(predictors_a, covars2), collapse = " + "), sep = " ~ ")), weights = weights, clusters = location_town)) %>% pool %>% summary(conf.int = TRUE) %>% mutate(outcome = "rtrn_head_ever")
reg2b = with(df, lm_robust(as.formula(paste("rtrn_head_ever", paste(c(predictors_b, covars2), collapse = " + "), sep = " ~ ")), weights = weights, clusters = location_town)) %>% pool %>% summary(conf.int = TRUE) %>% mutate(outcome = "rtrn_head_ever")
reg2 = reg2b %>% bind_rows(reg2a %>% filter(term == "index_safety_syria"))

reg3a = with(df, lm_robust(as.formula(paste("prep_pca", paste(c(predictors_a, covars2), collapse = " + "), sep = " ~ ")), weights = weights, clusters = location_town)) %>% pool %>% summary(conf.int = TRUE) %>% mutate(outcome = "prep_pca")
reg3b = with(df, lm_robust(as.formula(paste("prep_pca", paste(c(predictors_b, covars2), collapse = " + "), sep = " ~ ")), weights = weights, clusters = location_town)) %>% pool %>% summary(conf.int = TRUE) %>% mutate(outcome = "prep_pca")
reg3 = reg3b %>% bind_rows(reg3a %>% filter(term == "index_safety_syria"))

reg4a = with(df, lm_robust(as.formula(paste("expect2yr_syria", paste(c(predictors_a, covars2), collapse = " + "), sep = " ~ ")), weights = weights, clusters = location_town)) %>% pool %>% summary(conf.int = TRUE) %>% mutate(outcome = "expect2yr_syria")
reg4b = with(df, lm_robust(as.formula(paste("expect2yr_syria", paste(c(predictors_b, covars2), collapse = " + "), sep = " ~ ")), weights = weights, clusters = location_town)) %>% pool %>% summary(conf.int = TRUE) %>% mutate(outcome = "expect2yr_syria")
reg4 = reg4b %>% bind_rows(reg4a %>% filter(term == "index_safety_syria"))

dfplot3 = reg1 %>% bind_rows(reg2) %>% bind_rows(reg3) %>% bind_rows(reg4)

dfplot3 = dfplot3 %>% filter(term %in% names(indices))

dfplot3$term %>% unique

dfplot3 = dfplot3 %>% mutate(bold = "plain") %>% 
  mutate(term = factor(term, levels = c("index_safety_syria", "index_regime_control",
                                        "index_econ_syria", "index_services_syria", "index_networks_syria",
                                        "index_econ_leb",
                                        "index_services_leb", "index_networks_lebanon", "index_soc_leb",
                                        "index_legal_leb", "index_mobility", "index_family",
                                        "index_info_quality"))) %>% 
  mutate(term = recode(term, "index_econ_leb" = "Economic well-being (Leb.)",
                       "index_soc_leb" = "Social well-being (Leb.)",
                       "index_services_leb" = "Services (Leb.)",
                       "index_legal_leb" = "Legal conditions (Leb.)",
                       "index_safety_syria" = "Safety (Syr.)",
                       "index_econ_syria" = "Economic well-being (Syr.)",
                       "index_services_syria" = "Services (Syr.)",
                       "index_networks_syria" = "Networks (Syr.)",
                       "index_mobility" = "Log travel distance", 
                       "index_family" = "Log household size",
                       "index_networks_lebanon" = "Networks (Leb.)",
                       "index_info_quality" = "Confidence in information",
                       "index_regime_control" = "Regime control (Syr.)")) %>% 
  #mutate(term = fct_rev(term)) %>% 
  mutate(outcome = recode(outcome, "rtrn_head_12" = "Return in 12 months",
                          "rtrn_head_ever" = "Return ever",
                          "prep_pca" = "Prepare to return", 
                          "expect2yr_syria" = "Return in 2 years"),
         outcome = factor(outcome, levels = c("Return in 12 months", "Prepare to return", "Return in 2 years", 
                                              "Return ever"))
  )


reg1 = dfplot3 %>% filter(outcome == "Return in 12 months")
reg2 = dfplot3 %>% filter(outcome == "Prepare to return")

mod5 = createTexreg(coef.names = as.character(reg1$term), coef = reg1$estimate, se = reg1$std.error, pvalue = reg1$p.value, model.name = "12 months (Services)")
mod6 = createTexreg(coef.names = as.character(reg2$term), coef = reg2$estimate, se = reg2$std.error, pvalue = reg2$p.value, model.name = "Prepare (Services)")


# Excluding locality fixed effects -----------------------------------------------------------
predictors = names(indices)[2:length(indices)]

covars_nofe = covars[covars!="location_district"]
covars_nofe2 = covars2[covars2!="location_district"]

for(i in 1:length(predictors)){
  if(predictors[i] != "index_mobility"){
    assign(paste0("reg", i), with(df, lm_robust(as.formula(paste("rtrn_head_12", paste(c(predictors[i], covars_nofe), collapse = " + "), sep = " ~ ")), clusters = location_town, weights = weights)) %>% pool() %>% summary(conf.int = TRUE) %>% mutate(outcome = "rtrn_head_12"))
  }else{
    assign(paste0("reg", i), with(df, lm_robust(as.formula(paste("rtrn_head_12", paste(c(predictors[i], covars_nofe2), collapse = " + "), sep = " ~ ")), clusters = location_town, weights = weights)) %>% pool() %>% summary(conf.int = TRUE) %>% mutate(outcome = "rtrn_head_12"))
  }
}
df$syr_origin2 %>% class
df_plot = reg1

for(i in 2:length(predictors)){
  df_plot = df_plot %>% bind_rows(eval(parse(text = paste0("reg", i))))
}


for(i in 1:length(predictors)){
  if(predictors[i] != "index_mobility"){
    assign(paste0("reg", i), with(df, lm_robust(as.formula(paste("rtrn_head_ever", paste(c(predictors[i], covars_nofe), collapse = " + "), sep = " ~ ")), clusters = location_town, weights = weights)) %>% pool() %>% summary(conf.int = TRUE) %>% mutate(outcome = "rtrn_head_ever"))
  }else{
    assign(paste0("reg", i), with(df, lm_robust(as.formula(paste("rtrn_head_ever", paste(c(predictors[i], covars_nofe2), collapse = " + "), sep = " ~ ")), clusters = location_town, weights = weights)) %>% pool() %>% summary(conf.int = TRUE) %>% mutate(outcome = "rtrn_head_ever"))
  }
}

for(i in 1:length(predictors)){
  df_plot = df_plot %>% bind_rows(eval(parse(text = paste0("reg", i))))
}



for(i in 1:length(predictors)){
  if(predictors[i] != "index_mobility"){
    assign(paste0("reg", i), with(df, lm_robust(as.formula(paste("prep_pca", paste(c(predictors[i], covars_nofe), collapse = " + "), sep = " ~ ")), clusters = location_town, weights = weights)) %>% pool() %>% summary(conf.int = TRUE) %>% mutate(outcome = "prep_pca"))
  }else{
    assign(paste0("reg", i), with(df, lm_robust(as.formula(paste("prep_pca", paste(c(predictors[i], covars_nofe2), collapse = " + "), sep = " ~ ")), clusters = location_town, weights = weights)) %>% pool() %>% summary(conf.int = TRUE) %>% mutate(outcome = "prep_pca"))
  }
}

for(i in 1:length(predictors)){
  df_plot = df_plot %>% bind_rows(eval(parse(text = paste0("reg", i))))
}


for(i in 1:length(predictors)){
  if(predictors[i] != "index_mobility"){
    assign(paste0("reg", i), with(df, lm_robust(as.formula(paste("expect2yr_syria", paste(c(predictors[i], covars_nofe), collapse = " + "), sep = " ~ ")), clusters = location_town, weights = weights)) %>% pool() %>% summary(conf.int = TRUE) %>% mutate(outcome = "expect2yr_syria"))
  }else{
    assign(paste0("reg", i), with(df, lm_robust(as.formula(paste("expect2yr_syria", paste(c(predictors[i], covars_nofe2), collapse = " + "), sep = " ~ ")), clusters = location_town, weights = weights)) %>% pool() %>% summary(conf.int = TRUE) %>% mutate(outcome = "expect2yr_syria"))
  }
}

for(i in 1:length(predictors)){
  df_plot = df_plot %>% bind_rows(eval(parse(text = paste0("reg", i))))
}

df_plot = df_plot %>% filter(term %in% names(indices))

df_plot$term %>% unique

df_plot2 = df_plot %>% mutate(bold = "plain") %>% 
  mutate(term = factor(term, levels = c("index_safety_syria", "index_regime_control",
                                        "index_econ_syria", "index_services_syria", "index_networks_syria",
                                        "index_econ_leb",
                                        "index_services_leb", "index_networks_lebanon", "index_soc_leb",
                                        "index_legal_leb", "index_mobility", "index_family",
                                        "index_info_quality"))) %>% 
  mutate(term = recode(term, "index_econ_leb" = "Economic well-being (Leb.)",
                       "index_soc_leb" = "Social well-being (Leb.)",
                       "index_services_leb" = "Services (Leb.)",
                       "index_legal_leb" = "Legal conditions (Leb.)",
                       "index_safety_syria" = "Safety (Syr.)",
                       "index_econ_syria" = "Economic well-being (Syr.)",
                       "index_services_syria" = "Services (Syr.)",
                       "index_networks_syria" = "Networks (Syr.)",
                       "index_mobility" = "Log travel distance", 
                       "index_family" = "Log household size",
                       "index_networks_lebanon" = "Networks (Leb.)",
                       "index_info_quality" = "Confidence in information",
                       "index_regime_control" = "Regime control (Syr.)")) %>% 
  #mutate(term = fct_rev(term)) %>% 
  mutate(outcome = recode(outcome, "rtrn_head_12" = "Return in 12 months",
                          "rtrn_head_ever" = "Return ever",
                          "prep_pca" = "Prepare to return", 
                          "expect2yr_syria" = "Return in 2 years"),
         outcome = factor(outcome, levels = c("Return in 12 months", "Prepare to return", "Return in 2 years", 
                                              "Return ever"))
  )


df_plot2 = df_plot2 %>% arrange(outcome, term)
df_plot2 %>% filter(outcome == "Prepare to return") %>% select(term, estimate)
df_plot2 = df_plot2 %>% filter(!is.na(term))


reg1 = df_plot2 %>% filter(outcome == "Return in 12 months")
reg2 = df_plot2 %>% filter(outcome == "Prepare to return")

mod7 = createTexreg(coef.names = as.character(reg1$term), coef = reg1$estimate, se = reg1$std.error, pvalue = reg1$p.value, model.name = "12 months (No FEs)")
mod8 = createTexreg(coef.names = as.character(reg2$term), coef = reg2$estimate, se = reg2$std.error, pvalue = reg2$p.value, model.name = "Prepare (No FEs)")



# All indices appendix ----------------------------------------------------
predictors = names(indices)[2:length(indices)]

  reg1 = with(df, lm_robust(as.formula(paste("rtrn_head_12", paste(c(predictors, covars2), collapse = " + "), sep = " ~ ")), weights = weights, clusters = location_town)) %>% pool %>% summary(conf.int = TRUE) %>% mutate(outcome = "rtrn_head_12")
reg2 = with(df, lm_robust(as.formula(paste("rtrn_head_ever", paste(c(predictors, covars2), collapse = " + "), sep = " ~ ")), weights = weights, clusters = location_town)) %>% pool %>% summary(conf.int = TRUE) %>% mutate(outcome = "rtrn_head_ever")
reg3 = with(df, lm_robust(as.formula(paste("prep_pca", paste(c(predictors, covars2), collapse = " + "), sep = " ~ ")), weights = weights, clusters = location_town)) %>% pool %>% summary(conf.int = TRUE) %>% mutate(outcome = "prep_pca")
reg4 = with(df, lm_robust(as.formula(paste("expect2yr_syria", paste(c(predictors, covars2), collapse = " + "), sep = " ~ ")), weights = weights, clusters = location_town)) %>% pool %>% summary(conf.int = TRUE) %>% mutate(outcome = "expect2yr_syria")

dfplot2 = reg1 %>% bind_rows(reg2) %>% bind_rows(reg3) %>% bind_rows(reg4)

dfplot2 = dfplot2 %>% filter(term %in% names(indices))

dfplot2 = dfplot2 %>% mutate(bold = "plain") %>% 
  mutate(term = factor(term, levels = c("index_safety_syria", "index_regime_control",
                                        "index_econ_syria", "index_services_syria", "index_networks_syria",
                                        "index_econ_leb",
                                        "index_services_leb", "index_networks_lebanon", "index_soc_leb",
                                        "index_legal_leb", "index_mobility", "index_family",
                                        "index_info_quality"))) %>% 
  mutate(term = recode(term, "index_econ_leb" = "Economic well-being (Leb.)",
                       "index_soc_leb" = "Social well-being (Leb.)",
                       "index_services_leb" = "Services (Leb.)",
                       "index_legal_leb" = "Legal conditions (Leb.)",
                       "index_safety_syria" = "Safety (Syr.)",
                       "index_econ_syria" = "Economic well-being (Syr.)",
                       "index_services_syria" = "Services (Syr.)",
                       "index_networks_syria" = "Networks (Syr.)",
                       "index_mobility" = "Log travel distance", 
                       "index_family" = "Log household size",
                       "index_networks_lebanon" = "Networks (Leb.) ",
                       "index_info_quality" = "Confidence in information",
                       "index_regime_control" = "Regime control (Syr.)")) %>% 
  #mutate(term = fct_rev(term)) %>% 
  mutate(outcome = recode(outcome, "rtrn_head_12" = "Return in 12 months",
                          "rtrn_head_ever" = "Return ever",
                          "prep_pca" = "Prepare to return", 
                          "expect2yr_syria" = "Return in 2 years"),
         outcome = factor(outcome, levels = c("Return in 12 months", "Prepare to return", "Return in 2 years", 
                                              "Return ever"))
  )


dfplot2 = dfplot2 %>% arrange(outcome, term)


reg1 = dfplot2 %>% filter(outcome == "Return in 12 months")
reg2 = dfplot2 %>% filter(outcome == "Prepare to return")

mod9 = createTexreg(coef.names = as.character(reg1$term), coef = reg1$estimate, se = reg1$std.error, pvalue = reg1$p.value, model.name = "12 months (Safety + Services)")
mod10 = createTexreg(coef.names = as.character(reg2$term), coef = reg2$estimate, se = reg2$std.error, pvalue = reg2$p.value, model.name = "Prepare (Safety + Services)")








# Household return intentions ---------------------------------------------
df_long$outcome_hh = ifelse(df_long$rtrn_hh_12mon == 1 | df_long$rtrn_12mon == 1, 1, 0)
df_long$outcome_hh %>% table
df_long$rtrn_12mon %>% table
df_long$rtrn_hh_12mon %>% table

df_long$outcome_head2 = ifelse(df_long$rtrn_12mon %in% c(1, 2), 1, 0)

df = as.mids(df_long)

indices = df$data %>% 
  select(c(response_num, names(df$data)[str_detect(names(df$data), "^index")]))

covars = c("syr_origin", "syr_origin_urban", "location_its", "sick", "head_educ2", "toddler", "elderly", "female_headed_hh", "location_district", "hezb_control")
covars2 = c("syr_origin_urban", "location_its", "sick", "head_educ2", "toddler", "elderly", "female_headed_hh", "hezb_control")


#All indices
predictors = names(indices)[2:length(indices)]
predictors_a = predictors[which(predictors!="index_services_syria")]
predictors_b = predictors[which(predictors!="index_safety_syria")]

reg1a = with(df, lm_robust(as.formula(paste("outcome_hh", paste(c(predictors_a, covars2), collapse = " + "), sep = " ~ ")), weights = weights, clusters = location_town)) %>% pool %>% summary(conf.int = TRUE) %>% mutate(outcome = "outcome_hh")
reg1b = with(df, lm_robust(as.formula(paste("outcome_hh", paste(c(predictors_b, covars2), collapse = " + "), sep = " ~ ")), weights = weights, clusters = location_town)) %>% pool %>% summary(conf.int = TRUE) %>% mutate(outcome = "outcome_hh")
reg1 = reg1a %>% bind_rows(reg1b %>% filter(term == "index_services_syria"))

reg2a = with(df, lm_robust(as.formula(paste("outcome_head2", paste(c(predictors_a, covars2), collapse = " + "), sep = " ~ ")), weights = weights, clusters = location_town)) %>% pool %>% summary(conf.int = TRUE) %>% mutate(outcome = "outcome_head2")
reg2b = with(df, lm_robust(as.formula(paste("outcome_head2", paste(c(predictors_b, covars2), collapse = " + "), sep = " ~ ")), weights = weights, clusters = location_town)) %>% pool %>% summary(conf.int = TRUE) %>% mutate(outcome = "outcome_head2")
reg2 = reg2a %>% bind_rows(reg2b %>% filter(term == "index_services_syria"))

dfplot2 = reg1 %>% bind_rows(reg2)

dfplot2 = dfplot2 %>% filter(term %in% names(indices))

dfplot2 = dfplot2 %>% mutate(bold = "plain") %>% 
  mutate(term = factor(term, levels = c("index_safety_syria", "index_regime_control",
                                        "index_econ_syria", "index_services_syria", "index_networks_syria",
                                        "index_econ_leb",
                                        "index_services_leb", "index_networks_lebanon", "index_soc_leb",
                                        "index_legal_leb", "index_mobility", "index_family",
                                        "index_info_quality"))) %>% 
  mutate(term = recode(term, "index_econ_leb" = "Economic well-being (Leb.)",
                       "index_soc_leb" = "Social well-being (Leb.)",
                       "index_services_leb" = "Services (Leb.)",
                       "index_legal_leb" = "Legal conditions (Leb.)",
                       "index_safety_syria" = "Safety (Syr.)",
                       "index_econ_syria" = "Economic well-being (Syr.)",
                       "index_services_syria" = "Services (Syr.)",
                       "index_networks_syria" = "Networks (Syr.)",
                       "index_mobility" = "Log travel distance", 
                       "index_family" = "Log household size",
                       "index_networks_lebanon" = "Networks (Leb.) ",
                       "index_info_quality" = "Confidence in information",
                       "index_regime_control" = "Regime control (Syr.)")) %>% 
  #mutate(term = fct_rev(term)) %>% 
  mutate(outcome = recode(outcome, "outcome_hh" = "Household member return",
                          "outcome_head2" = "Head return (counting uncertain)"),
         outcome = factor(outcome, levels = c("Household member return", "Head return (counting uncertain)")))


reg1 = dfplot2 %>% filter(outcome == "Household member return")
reg2 = dfplot2 %>% filter(outcome == "Head return (counting uncertain)")

mod11 = createTexreg(coef.names = as.character(reg1$term), coef = reg1$estimate, se = reg1$std.error, pvalue = reg1$p.value, model.name = "Household member return")
mod12 = createTexreg(coef.names = as.character(reg2$term), coef = reg2$estimate, se = reg2$std.error, pvalue = reg2$p.value, model.name = "Head of HH return (counting uncertain)")



# Additive indices --------------------------------------------------------
## Individual indices --------------------------------------------------------
indices = df$data %>% 
  select(c(response_num, names(df$data)[str_detect(names(df$data), "^additive")]))

predictors = names(indices)[2:length(indices)]

for(i in 1:length(predictors)){
  if(predictors[i] != "additive_mobility"){
    assign(paste0("reg", i), with(df, lm_robust(as.formula(paste("rtrn_head_12", paste(c(predictors[i], covars), collapse = " + "), sep = " ~ ")), clusters = location_town, weights = weights)) %>% pool() %>% summary(conf.int = TRUE) %>% mutate(outcome = "rtrn_head_12"))
  }else{
    assign(paste0("reg", i), with(df, lm_robust(as.formula(paste("rtrn_head_12", paste(c(predictors[i], covars2), collapse = " + "), sep = " ~ ")), clusters = location_town, weights = weights)) %>% pool() %>% summary(conf.int = TRUE) %>% mutate(outcome = "rtrn_head_12"))
  }
}

dfplot1 = reg1

for(i in 2:length(predictors)){
  dfplot1 = dfplot1 %>% bind_rows(eval(parse(text = paste0("reg", i))))
}

for(i in 1:length(predictors)){
  if(predictors[i] != "additive_mobility"){
    assign(paste0("reg", i), with(df, lm_robust(as.formula(paste("rtrn_head_ever", paste(c(predictors[i], covars), collapse = " + "), sep = " ~ ")), clusters = location_town, weights = weights)) %>% pool() %>% summary(conf.int = TRUE) %>% mutate(outcome = "rtrn_head_ever"))
  }else{
    assign(paste0("reg", i), with(df, lm_robust(as.formula(paste("rtrn_head_ever", paste(c(predictors[i], covars2), collapse = " + "), sep = " ~ ")), clusters = location_town, weights = weights)) %>% pool() %>% summary(conf.int = TRUE) %>% mutate(outcome = "rtrn_head_ever"))
  }
}

for(i in 1:length(predictors)){
  dfplot1 = dfplot1 %>% bind_rows(eval(parse(text = paste0("reg", i))))
}


for(i in 1:length(predictors)){
  if(predictors[i] != "additive_mobility"){
    assign(paste0("reg", i), with(df, lm_robust(as.formula(paste("prep_pca", paste(c(predictors[i], covars), collapse = " + "), sep = " ~ ")), clusters = location_town, weights = weights)) %>% pool() %>% summary(conf.int = TRUE) %>% mutate(outcome = "prep_pca"))
  }else{
    assign(paste0("reg", i), with(df, lm_robust(as.formula(paste("prep_pca", paste(c(predictors[i], covars2), collapse = " + "), sep = " ~ ")), clusters = location_town, weights = weights)) %>% pool() %>% summary(conf.int = TRUE) %>% mutate(outcome = "prep_pca"))
  }
}

for(i in 1:length(predictors)){
  dfplot1 = dfplot1 %>% bind_rows(eval(parse(text = paste0("reg", i))))
}


for(i in 1:length(predictors)){
  if(predictors[i] != "additive_mobility"){
    assign(paste0("reg", i), with(df, lm_robust(as.formula(paste("expect2yr_syria", paste(c(predictors[i], covars), collapse = " + "), sep = " ~ ")), clusters = location_town, weights = weights)) %>% pool() %>% summary(conf.int = TRUE) %>% mutate(outcome = "expect2yr_syria"))
  }else{
    assign(paste0("reg", i), with(df, lm_robust(as.formula(paste("expect2yr_syria", paste(c(predictors[i], covars2), collapse = " + "), sep = " ~ ")), clusters = location_town, weights = weights)) %>% pool() %>% summary(conf.int = TRUE) %>% mutate(outcome = "expect2yr_syria"))
  }
}

for(i in 1:length(predictors)){
  dfplot1 = dfplot1 %>% bind_rows(eval(parse(text = paste0("reg", i))))
}

dfplot1 = dfplot1 %>% filter(term %in% names(indices))

dfplot1 = dfplot1 %>% mutate(bold = "plain") %>% 
  mutate(term = factor(term, levels = c("additive_safety_syria", 
                                        "additive_econ_syria", "additive_services_syria", "additive_networks_syria",
                                        "additive_econ_leb",
                                        "additive_services_leb", "additive_networks_lebanon", "additive_soc_leb",
                                        "additive_legal_leb", "additive_mobility", "additive_family",
                                        "additive_info_quality"))) %>% 
  mutate(term = recode(term, "additive_econ_leb" = "Economic well-being (Leb.)",
                       "additive_soc_leb" = "Social well-being (Leb.)",
                       "additive_services_leb" = "Services (Leb.)",
                       "additive_legal_leb" = "Legal conditions (Leb.)",
                       "additive_safety_syria" = "Safety (Syr.)",
                       "additive_econ_syria" = "Economic well-being (Syr.)",
                       "additive_services_syria" = "Services (Syr.)",
                       "additive_networks_syria" = "Networks (Syr.)",
                       "additive_mobility" = "Log travel distance", 
                       "additive_family" = "Log household size",
                       "additive_networks_lebanon" = "Networks (Leb.) ",
                       "additive_info_quality" = "Confidence in information")) %>% 
  #mutate(term = fct_rev(term)) %>% 
  mutate(outcome = recode(outcome, "rtrn_head_12" = "Return in 12 months",
                          "rtrn_head_ever" = "Return ever",
                          "prep_pca" = "Prepare to return", 
                          "expect2yr_syria" = "Return in 2 years"),
         outcome = factor(outcome, levels = c("Return in 12 months", "Prepare to return", "Return in 2 years", 
                                              "Return ever"))
  )

dfplot1 = dfplot1 %>% arrange(outcome, term)
dfplot1 = dfplot1 %>% filter(!is.na(term))

reg1 = dfplot1 %>% filter(outcome == "Return in 12 months")
reg2 = dfplot1 %>% filter(outcome == "Prepare to return")

mod13 = createTexreg(coef.names = as.character(reg1$term), coef = reg1$estimate, se = reg1$std.error, pvalue = reg1$p.value, model.name = "12 months (Ind. Indices)")
mod14 = createTexreg(coef.names = as.character(reg2$term), coef = reg2$estimate, se = reg2$std.error, pvalue = reg2$p.value, model.name = "Prepare (Ind. Indices)")


## Regressions all indices-------------------------------------------------------------
predictors = names(indices)[2:length(indices)]
predictors_a = predictors[which(predictors!="additive_services_syria")]
predictors_b = predictors[which(predictors!="additive_safety_syria")]

reg1a = with(df, lm_robust(as.formula(paste("rtrn_head_12", paste(c(predictors_a, covars2), collapse = " + "), sep = " ~ ")), weights = weights, clusters = location_town)) %>% pool %>% summary(conf.int = TRUE) %>% mutate(outcome = "rtrn_head_12")
reg1b = with(df, lm_robust(as.formula(paste("rtrn_head_12", paste(c(predictors_b, covars2), collapse = " + "), sep = " ~ ")), weights = weights, clusters = location_town)) %>% pool %>% summary(conf.int = TRUE) %>% mutate(outcome = "rtrn_head_12")
reg1 = reg1a %>% bind_rows(reg1b %>% filter(term == "additive_services_syria"))

reg2a = with(df, lm_robust(as.formula(paste("rtrn_head_ever", paste(c(predictors_a, covars2), collapse = " + "), sep = " ~ ")), weights = weights, clusters = location_town)) %>% pool %>% summary(conf.int = TRUE) %>% mutate(outcome = "rtrn_head_ever")
reg2b = with(df, lm_robust(as.formula(paste("rtrn_head_ever", paste(c(predictors_b, covars2), collapse = " + "), sep = " ~ ")), weights = weights, clusters = location_town)) %>% pool %>% summary(conf.int = TRUE) %>% mutate(outcome = "rtrn_head_ever")
reg2 = reg2a %>% bind_rows(reg2b %>% filter(term == "additive_services_syria"))

reg3a = with(df, lm_robust(as.formula(paste("prep_pca", paste(c(predictors_a, covars2), collapse = " + "), sep = " ~ ")), weights = weights, clusters = location_town)) %>% pool %>% summary(conf.int = TRUE) %>% mutate(outcome = "prep_pca")
reg3b = with(df, lm_robust(as.formula(paste("prep_pca", paste(c(predictors_b, covars2), collapse = " + "), sep = " ~ ")), weights = weights, clusters = location_town)) %>% pool %>% summary(conf.int = TRUE) %>% mutate(outcome = "prep_pca")
reg3 = reg3a %>% bind_rows(reg3b %>% filter(term == "additive_services_syria"))

reg4a = with(df, lm_robust(as.formula(paste("expect2yr_syria", paste(c(predictors_a, covars2), collapse = " + "), sep = " ~ ")), weights = weights, clusters = location_town)) %>% pool %>% summary(conf.int = TRUE) %>% mutate(outcome = "expect2yr_syria")
reg4b = with(df, lm_robust(as.formula(paste("expect2yr_syria", paste(c(predictors_b, covars2), collapse = " + "), sep = " ~ ")), weights = weights, clusters = location_town)) %>% pool %>% summary(conf.int = TRUE) %>% mutate(outcome = "expect2yr_syria")
reg4 = reg4a %>% bind_rows(reg4b %>% filter(term == "additive_services_syria"))

dfplot2 = reg1 %>% bind_rows(reg2) %>% bind_rows(reg3) %>% bind_rows(reg4)

dfplot2 = dfplot2 %>% filter(term %in% names(indices))

dfplot2 = dfplot2 %>% mutate(bold = "plain") %>% 
  mutate(term = factor(term, levels = c("additive_safety_syria",
                                        "additive_econ_syria", "additive_services_syria", "additive_networks_syria",
                                        "additive_econ_leb",
                                        "additive_services_leb", "additive_networks_lebanon", "additive_soc_leb",
                                        "additive_legal_leb", "additive_mobility", "additive_family",
                                        "additive_info_quality"))) %>% 
  mutate(term = recode(term, "additive_econ_leb" = "Economic well-being (Leb.)",
                       "additive_soc_leb" = "Social well-being (Leb.)",
                       "additive_services_leb" = "Services (Leb.)",
                       "additive_legal_leb" = "Legal conditions (Leb.)",
                       "additive_safety_syria" = "Safety (Syr.)",
                       "additive_econ_syria" = "Economic well-being (Syr.)",
                       "additive_services_syria" = "Services (Syr.)",
                       "additive_networks_syria" = "Networks (Syr.)",
                       "additive_mobility" = "Log travel distance", 
                       "additive_family" = "Log household size",
                       "additive_networks_lebanon" = "Networks (Leb.) ",
                       "additive_info_quality" = "Confidence in information")) %>% 
  #mutate(term = fct_rev(term)) %>% 
  mutate(outcome = recode(outcome, "rtrn_head_12" = "Return in 12 months",
                          "rtrn_head_ever" = "Return ever",
                          "prep_pca" = "Prepare to return", 
                          "expect2yr_syria" = "Return in 2 years"),
         outcome = factor(outcome, levels = c("Return in 12 months", "Prepare to return", "Return in 2 years", 
                                              "Return ever"))
  )


reg1 = dfplot2 %>% filter(outcome == "Return in 12 months")
reg2 = dfplot2 %>% filter(outcome == "Prepare to return")

mod15 = createTexreg(coef.names = as.character(reg1$term), coef = reg1$estimate, se = reg1$std.error, pvalue = reg1$p.value, model.name = "12 months (All Indices)")
mod16 = createTexreg(coef.names = as.character(reg2$term), coef = reg2$estimate, se = reg2$std.error, pvalue = reg2$p.value, model.name = "Prepare (All Indices)")


# Individual regressions with additional covariates -----------------------------------------------------------
df_long$need_school = 1 - df_long$no_need_school
df_long$curfews_now = 1 - df_long$no_curfews_now
df = as.mids(df_long)

reg1 = with(df, lm_robust(as.formula(paste("rtrn_head_12", paste(c("index_econ_leb", covars, "length_stay", "curfews_now", "need_school", "leban_relatives2"), collapse = "+"), sep = " ~ ")), clusters = location_town, weights = weights)) %>% pool() %>% summary(conf.int = TRUE) %>% mutate(outcome = "rtrn_head_12")
reg2 = with(df, lm_robust(as.formula(paste("rtrn_head_12", paste(c("index_soc_leb", covars, "hh_income", "hh_inc_source_aid", "need_school", "leban_relatives2"), collapse = "+"), sep = " ~ ")), clusters = location_town, weights = weights)) %>% pool() %>% summary(conf.int = TRUE) %>% mutate(outcome = "rtrn_head_12")
reg3 = with(df, lm_robust(as.formula(paste("rtrn_head_12", paste(c("index_services_leb", covars, "hh_income", "hh_inc_source_aid", "length_stay", "curfews_now", "leban_relatives2"), collapse = "+"), sep = " ~ ")), clusters = location_town, weights = weights)) %>% pool() %>% summary(conf.int = TRUE) %>% mutate(outcome = "rtrn_head_12")
reg4 = with(df, lm_robust(as.formula(paste("rtrn_head_12", paste(c("index_legal_leb", covars, "hh_income", "hh_inc_source_aid", "length_stay", "curfews_now", "need_school", "leban_relatives2"), collapse = "+"), sep = " ~ ")), clusters = location_town, weights = weights)) %>% pool() %>% summary(conf.int = TRUE) %>% mutate(outcome = "rtrn_head_12")
reg5 = with(df, lm_robust(as.formula(paste("rtrn_head_12", paste(c("index_networks_lebanon", covars, "hh_income", "hh_inc_source_aid", "length_stay", "curfews_now", "need_school"), collapse = "+"), sep = " ~ ")), clusters = location_town, weights = weights)) %>% pool() %>% summary(conf.int = TRUE) %>% mutate(outcome = "rtrn_head_12")
reg6 = with(df, lm_robust(as.formula(paste("rtrn_head_12", paste(c("index_safety_syria", covars, "hh_income", "hh_inc_source_aid", "length_stay", "curfews_now", "need_school", "leban_relatives2"), collapse = "+"), sep = " ~ ")), clusters = location_town, weights = weights)) %>% pool() %>% summary(conf.int = TRUE) %>% mutate(outcome = "rtrn_head_12")
reg7 = with(df, lm_robust(as.formula(paste("rtrn_head_12", paste(c("index_regime_control", covars, "hh_income", "hh_inc_source_aid", "length_stay", "curfews_now", "need_school", "leban_relatives2"), collapse = "+"), sep = " ~ ")), clusters = location_town, weights = weights)) %>% pool() %>% summary(conf.int = TRUE) %>% mutate(outcome = "rtrn_head_12")
reg8 = with(df, lm_robust(as.formula(paste("rtrn_head_12", paste(c("index_econ_syria", covars, "hh_income", "hh_inc_source_aid", "length_stay", "curfews_now", "need_school", "leban_relatives2"), collapse = "+"), sep = " ~ ")), clusters = location_town, weights = weights)) %>% pool() %>% summary(conf.int = TRUE) %>% mutate(outcome = "rtrn_head_12")
reg9 = with(df, lm_robust(as.formula(paste("rtrn_head_12", paste(c("index_services_syria", covars, "hh_income", "hh_inc_source_aid", "length_stay", "curfews_now", "need_school", "leban_relatives2"), collapse = "+"), sep = " ~ ")), clusters = location_town, weights = weights)) %>% pool() %>% summary(conf.int = TRUE) %>% mutate(outcome = "rtrn_head_12")
reg10 = with(df, lm_robust(as.formula(paste("rtrn_head_12", paste(c("index_networks_syria", covars, "hh_income", "hh_inc_source_aid", "length_stay", "curfews_now", "need_school", "leban_relatives2"), collapse = "+"), sep = " ~ ")), clusters = location_town, weights = weights)) %>% pool() %>% summary(conf.int = TRUE) %>% mutate(outcome = "rtrn_head_12")
reg11 = with(df, lm_robust(as.formula(paste("rtrn_head_12", paste(c("index_info_quality", covars, "hh_income", "hh_inc_source_aid", "length_stay", "curfews_now", "need_school", "leban_relatives2"), collapse = "+"), sep = " ~ ")), clusters = location_town, weights = weights)) %>% pool() %>% summary(conf.int = TRUE) %>% mutate(outcome = "rtrn_head_12")
reg12 = with(df, lm_robust(as.formula(paste("rtrn_head_12", paste(c("index_mobility", covars2, "hh_income", "hh_inc_source_aid", "length_stay", "curfews_now", "need_school", "leban_relatives2"), collapse = "+"), sep = " ~ ")), clusters = location_town, weights = weights)) %>% pool() %>% summary(conf.int = TRUE) %>% mutate(outcome = "rtrn_head_12")
reg13 = with(df, lm_robust(as.formula(paste("rtrn_head_12", paste(c("index_family", covars, "hh_income", "hh_inc_source_aid", "length_stay", "curfews_now", "need_school", "leban_relatives2"), collapse = "+"), sep = " ~ ")), clusters = location_town, weights = weights)) %>% pool() %>% summary(conf.int = TRUE) %>% mutate(outcome = "rtrn_head_12")

dfplot1 = reg1 %>% bind_rows(reg2) %>% bind_rows(reg3) %>% bind_rows(reg4) %>% bind_rows(reg5) %>% bind_rows(reg6) %>% bind_rows(reg7) %>% bind_rows(reg8) %>% bind_rows(reg9) %>% bind_rows(reg10) %>% bind_rows(reg11) %>% bind_rows(reg12) %>% bind_rows(reg13)

dfplot1 = dfplot1 %>% mutate(bold = "plain") %>% 
  mutate(term = factor(term, levels = c("index_safety_syria", "index_regime_control",
                                        "index_econ_syria", "index_services_syria", "index_networks_syria",
                                        "index_econ_leb",
                                        "index_services_leb", "index_networks_lebanon", "index_soc_leb",
                                        "index_legal_leb", "index_mobility", "index_family",
                                        "index_info_quality"))) %>% 
  mutate(term = recode(term, "index_econ_leb" = "Economic well-being (Leb.)",
                       "index_soc_leb" = "Social well-being (Leb.)",
                       "index_services_leb" = "Services (Leb.)",
                       "index_legal_leb" = "Legal conditions (Leb.)",
                       "index_safety_syria" = "Safety (Syr.)",
                       "index_econ_syria" = "Economic well-being (Syr.)",
                       "index_services_syria" = "Services (Syr.)",
                       "index_networks_syria" = "Networks (Syr.)",
                       "index_mobility" = "Log travel distance", 
                       "index_family" = "Log household size",
                       "index_networks_lebanon" = "Networks (Leb.) ",
                       "index_info_quality" = "Confidence in information",
                       "index_regime_control" = "Regime control (Syr.)")) %>% 
  #mutate(term = fct_rev(term)) %>% 
  mutate(outcome = recode(outcome, "rtrn_head_12" = "Return in 12 months",
                          "rtrn_head_ever" = "Return ever",
                          "prep_pca" = "Prepare to return", 
                          "expect2yr_syria" = "Return in 2 years"),
         outcome = factor(outcome, levels = c("Return in 12 months", "Prepare to return", "Return in 2 years", 
                                              "Return ever"))
  )

dfplot1 = dfplot1 %>% arrange(outcome, term)
dfplot1 = dfplot1 %>% filter(!is.na(term))

mod17 = createTexreg(coef.names = as.character(dfplot1$term), coef = dfplot1$estimate, se = dfplot1$std.error, pvalue = dfplot1$p.value, model.name = "Return in 12 months")



# Export ------------------------------------------------------------------
texreg(list(mod1, mod2, mod3, mod4), file = "tables/robustness1.tex", stars = c(0.01, 0.05, 0.1), digits = 3, scalebox = .7, float.pos = "h!", caption = "Regression results using alternative outcomes. The first two models present the regression results for the return ever outcome 
       using one index per regression (model 1) and using all indices in the same regression (model 2). The last two models present the regression results for the return in 2 years outcome using one index per regression (model 3) and using all indices in the same regression (model 4).", label = "tab:robustness1", use.packages = F)

texreg(list(mod5, mod6, mod9, mod10, mod17), file = "tables/robustness2.tex", stars = c(0.01, 0.05, 0.1), digits = 3, scalebox = .7, float.pos = "h!", caption = "Robustness tests. The first two models include the services in Syria index in the all-indices regression (instead of the safety index)
       for the return in 12 months outcome (model 1) and the preparation to return outcome (model 2). The next two models present the analysis that includes running both the safety and services in Syria indices in the same regression. Note that there is a high correlation between the safety and services in Syria 
       indices. The third model uses the return in 12 months outcome and the fourth model uses the preparation to return outcome. The final model uses the main return in 12 months outcome with individual indices but uses additional covariates. Besides the ones mentioned in the paper, it controls for household income,
       an indicator for receiving aid, years in Lebanon, an indicator for a curfew in the locality, an indicator for age-school children, and an indicator for Lebanese relatives. Note that these additional controls are included only when the individual index in the regression does not contain any of these variables. 
       The list of variables in each index is included appendix section 3.", label = "tab:robustness2", use.packages = F)

texreg(list(mod7, mod8, mod11, mod12), file = "tables/robustness3.tex", stars = c(0.01, 0.05, 0.1), digits = 3, scalebox = .7, float.pos = "h!", caption = "Additional tests. The first two models re-run the main analysis without fixed effects for the return in 12 months outcome (model 1) and for the preparation to 
return outcome (model 2). The third models uses an alternative dependent variable: household member return intentions. This outcome is coded as 1 if any household member plans to return and 0 otherwise. The final model uses another outcome. The main return intentions in 12 month is coded as 1 if the respondent indicates their plan 
to return in 12 months and 0 otherwise. The last model uses the return in 12 months outcome but it is coded as 1 if the respondent indicates either that they plan to return in 12 months or they are uncertain about returning in 12 months and 0 otherwise.", label = "tab:robustness3", use.packages = F)

texreg(list(mod13, mod14, mod15, mod16), file = "tables/robustness4.tex", stars = c(0.01, 0.05, 0.1), digits = 3, scalebox = .7, float.pos = "h!", caption = "Additive indices. While the results in the paper use the first principal component of independent variables, here we construct the independent variables using additive indices. 
First, we standardize each variable to have zero mean and unit standard deviation. Then we summed these variables and standardized the sum again using the same method. The first two models present the results using one index in each regression for the return in 12 months outcome (model 1) and the preparation to return outcome (model 2). The last two 
models present the results using all the additive indices in one regression for the return in 12 months outcome (model 3) and the preparation to return outcome (model 4).", label = "tab:robustness4", use.packages = F)
